%                   SST(1^o)   SST(5^o)
% raw location      0          3
% up  hill          2          1

function [AT_subset,SST_subset] = SAT_SST_func_pair_CMIP6(lon,lat,P)

    dir_load     = SAT_SST_func_IO('CMIP6');
    file_save    = [dir_load,'CMIP6_tas_and_tos_1x1.mat'];  
    f            = matfile(file_save);
    
    if ismember(P.do_CMIP6_sens,[1 3])
        file_save    = [dir_load,'CMIP6_tos_5X5.mat'];  
        f2           = matfile(file_save);
    end
    
    N_mon        = 12;
    N_yr         = size(f.tas_hist(1,1,:,1),3) / N_mon;
    AT_subset    = nan(N_mon,N_yr,numel(lon));
    SST_subset   = nan(N_mon,N_yr,numel(lon));
    SST_subset_n = nan(N_mon,N_yr,numel(lon),4);
        
    % compute the new position of land stations ---------------------------
    if ismember(P.do_CMIP6_sens,[1 2])
        [lon2,lat2] = CDF_up_hill(lon,lat,1);
    end

    for ct_yr = 1:N_yr
        if rem(ct_yr+1849,10) == 0
            disp(['Pairing CMIP6 :: year ',num2str(ct_yr+1849)])
        end
        for ct_mon = 1:N_mon
            
            clear('tas_hist','tos_hist')
            id       = ct_mon + (ct_yr - 1) * 12;
            tas_hist = squeeze(f.tas_hist(:,:,id,:));
            if ismember(P.do_CMIP6_sens,[0 2])
                tos_hist = squeeze(f.tos_hist(:,:,id,:));     reso_sst = 1;
            elseif ismember(P.do_CMIP6_sens,[1 3])
                tos_hist = squeeze(f2.tos_hist(:,:,id,:));    reso_sst = 5;
            end

            for en = P.model_id

                if ismember(P.do_CMIP6_sens,[0 3])
                    AT_subset(ct_mon,ct_yr,:) = grd2pnt(lon,lat,tas_hist(:,:,en),1,1);
                elseif ismember(P.do_CMIP6_sens,[1 2])
                    AT_subset(ct_mon,ct_yr,:) = grd2pnt(lon2,lat2,tas_hist(:,:,en),1,1);
                end

                SST_subset(ct_mon,ct_yr,:) = ...
                       grd2pnt(lon,lat,tos_hist(:,:,en),reso_sst,reso_sst);

                SST_subset_n(ct_mon,ct_yr,:,1) = grd2pnt(...
                      lon+reso_sst,lat,tos_hist(:,:,en),reso_sst,reso_sst);
                SST_subset_n(ct_mon,ct_yr,:,2) = grd2pnt(...
                      lon-reso_sst,lat,tos_hist(:,:,en),reso_sst,reso_sst);
                SST_subset_n(ct_mon,ct_yr,:,3) = grd2pnt(...
                      lon,lat+reso_sst,tos_hist(:,:,en),reso_sst,reso_sst);
                SST_subset_n(ct_mon,ct_yr,:,4) = grd2pnt(...
                      lon,lat-reso_sst,tos_hist(:,:,en),reso_sst,reso_sst);
            end
        end
    end

    SST_subset_m                    = nanmean(SST_subset_n,4);
    SST_subset(isnan(SST_subset))   = SST_subset_m(isnan(SST_subset));
    
    AT_subset  = AT_subset(:,[P.yr_sub_st:P.yr_sub_ed]-1849,:);
    SST_subset = SST_subset(:,[P.yr_sub_st:P.yr_sub_ed]-1849,:);
    
end
